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In spite of the large number of papers appeared in the past which are devoted to the lattice 
Boltzmann (LB) methods, basic aspects of the theory still remain unchallenged. An unsolved theo- 
retical issue is related to the construction of a discrete kinetic theory which yields exactly the fluid 
equations, i.e., is non-asymptotic (here denoted as LB inverse kinetic theory). The purpose of this 
paper is theoretical and aims at developing an inverse kinetic approach of this type. In principle 
infinite solutions exist to this problem but the freedom can be exploited in order to meet important 
requirements. In particular, the discrete kinetic theory can be defined so that it yields exactly the 
fluid equation also for arbitrary non-equilibrium (but suitably smooth) kinetic distribution func- 
tions and arbitrarily close to the boundary of the fluid domain. This includes the specification 
of the kinetic initial and boundary conditions which are consistent with the initial and boundary 
conditions prescribed for the fluid fields. Other basic features are the arbitrariness of the "equi- 
librium" distribution function and the condition of positivity imposed on the kinetic distribution 
function. The latter can be achieved by imposing a suitable entropic principle, realized by means of 
a constant H-theorem. Unlike previous entropic LB methods the theorem can be obtained without 
functional constraints on the class of the initial distribution functions. As a basic consequence, the 
choice of the the entropy functional remains essentially arbitrary so that it can be identified with 
the Gibbs-Shannon entropy. Remarkably, this property is not affected by the particular choice of 
the kinetic equilibrium (to be assumed in all cases strictly positive). Hence, it applies also in the 
case of polynomial equilibria, usually adopted in customary LB approaches. We provide different 
possible realizations of the theory and asymptotic approximations which permit to determine the 
fluid equations with prescribed accuracy. As a result, asymptotic accuracy estimates of customary 
LB approaches and comparisons with the Chorin artificial compressibility method are discussed. 

PACS numbers: 47.27.Ak, 47.27.eb, 47.27.ed 
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1 - INTRODUCTION - INVERSE KINETIC 
THEORIES 



Basic issues concerning the foundations classical hy- 
drodynamics still remain unanswered. A remarkable as- 
pect is related the construction of inverse kinetic theo- 
ries (IKT) for hydrodynamic equations in which the fluid 
fields are identified with suitable moments of an appropri- 
ate kinetic probability distribution. The topic has been 
the subject of theoretical investigations both regarding 
the incompressible Navier-Stokes (NS) equations (INSE) 
[IB Sill SB] and the quantum hydrodynamic equations 
associated to the Schrodinger equation Q. The impor- 
tance of the IKT-approach for classical hydrodynamics 
goes beyond the academic interest. In fact, INSE rep- 
resent a mixture of hyperbolic and elliptic pde's, which 
are extremely hard to study both analytically and nu- 
merically. As such, their investigation represents a chal- 
lenge both for mathematical analysis and for computa- 
tional fluid dynamics. The discovery of IKT [l[ provides, 
however, a new starting point for the theoretical and nu- 
merical investigation of INSE. In fact, an inverse kinetic 



theory yields, by definition, an exact solver for the fluid 
equations: all the fluid fields, including the fluid pres- 
sure p(r, t) , are uniquely prescribed in terms of suitable 
momenta of the kinetic distribution function, solution 
of the kinetic equation. In the case of INSE this per- 
mits, in principle, to determine the evolution of the fluid 
fields without solving explicitly the Navier-Stokes equa- 
tion, nor the Poisson equations for the fluid pressure [|| . 
Previous IKT approaches @, d, 0, H, 0] have been based 
on continuous phase-space models. However, the inter- 
esting question arises whether similar concepts can be 
adopted also to the development of discrete inverse ki- 
netic theories based on the lattice Boltzmann (LB) the- 
ory. The goal of this investigation is to propose a novel 
LB theory for INSE, based on the development of an IKT 
with discrete velocities, here denoted as lattice Boltzmann 
inverse kinetic theory (LB-IKT). In this paper we intend 
to analyze the theoretical foundations and basic proper- 
ties of the new approach useful to display its relation- 
ship with previous CFD and lattice Boltzmann methods 
(LBM) for incompressible isothermal fluids. In particu- 
lar, we wish to prove that it delivers an inverse kinetic 
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theory, i.e., that it realizes an exact Navier-Stokes and 
Poisson solver. 



la - Motivations: difficulties with LBM's 

Despite the significant number of theoretical and nu- 
merical papers appeared in the literature in the last few 
years, the lattice Boltzmann method [1, 0, U, 12, 
ll3L 1 141 ] - among many others available in CFD - is prob- 
ably the one for which a complete understanding is not 
yet available. Although originated as an extension of 
the lattice gas automaton [13, 1 1 611 o r a special discrete 
form of the Boltzmann equation [l7j |. several aspects re- 
garding the very foundation of LB theory still remain to 
be clarified. Consequently, also the comparisons and ex- 
act relationship between the various lattice Boltzmann 
methods (LBM) and other CFD methods are made dif- 
ficult or, at least, not yet well understood. Needless to 
say, these comparisons are essential to assess the relative 
value (based on the characteristic computational com- 
plexity, accuracy and stability) of LBM and other CFD 
methods. In particular the relative performance of the 
numerical methods depend strongly on the characteris- 
tic spatial and time discretization scales, i.e., the minimal 
spatial and time scale lengths required by each numerical 
method to achieve a prescribed accuracy. On the other 
hand, most of the existing knowledge of the LBM's prop- 
erties originates from numerical benchmarks (see for ex- 
ample 18, 1^, 2~3])- Although these studies have demon- 
strated the LBM's accuracy in simulating fluid flows, few 
comparisons are available on the relative computational 
efficiency of the LBM and other CFD methods [HE]. 
The main reason [of these difficulties] is probably because 
current LBM's, rather than being exact Navier-Stokes 
solvers, are at most asymptotic ones (asymptotic LBM's), 
i.e., they depend on one or more infinitesimal parame- 
ters and recover INSE only in an approximate asymptotic 
sense. 

The motivations of this work are related to some of 
the basic features of customary LB theory representing, 
at the same time, assets and weaknesses. One of the 
main reasons of the popularity of the LB approach lays 
in its simplicity and in the fact that it provides an ap- 
proximate Poisson solver, i.e., it permits to advance in 
time the fluid fields without explicitly solving numeri- 
cally the Poisson equation for the fluid pressure. How- 
ever customary LB approaches can yield, at most, only 
asymptotic approximations for the fluid fields. This is 
because of two different reasons. The first one is the dif- 
ficulty in the precise definition of the kinetic boundary 
conditions in customary LBM's, since sufficiently close to 
the boundary the form of the distribution function pre- 
scribed by the boundary conditions is not generally con- 
sistent with hydrodynamic equations. The second reason 
is that the kinetic description adopted implies either the 



introduction of weak compressibility [1,0 [ll, 12, 13, 14\ 
or temperature [^ effects of the fluid or some sort of 
state equation for the fluid pressure [11]. These assump- 
tions, although physically plausible, appear unacceptable 
from the mathematical viewpoint since they represent a 
breaking of the exact fluid equations. 

Moreover, in the case of very small fluid viscosity 
customary LBM's may become inefficient as a conse- 
quence of the low-order approximations usually adopted 
and the possible presence of the numerical instabilities 
mentioned above. These accuracy limitations at low vis- 
cosities can usually be overcome only by imposing severe 
grid refinements and strong reductions of the size of the 
time step. This has the inevitable consequence of rais- 
ing significantly the level of computational complexity 
in customary LBM's (potentially much higher than that 
of so-called direct solution methods), which makes them 
inefficient or even potentially unsuitable for large-scale 
simulations in fluids. 

A fundamental issue is, therefore, related to the con- 
struction of more accurate, or higher-order, LBM's, ap- 
plicable for arbitrary values of the relevant physical 
(and asymptotic) parameters. However, the route which 
should permit to determine them is still uncertain, since 
the very existence of an underlying exact (and non- 
asymptotic) discrete kinetic theory-jjinalogous to the con- 



is not yet known. 
this should be 



;neory, 

tinuous inverse kinetic theor y I2L \l 
According to some authors [2J, [25, 
linked to the discretization of the Boltzmann equation, or 
to the possible introduction of weakly compressible and 
thermal flow models. However, the first approach is not 
only extremely hard to implement [2t| . since it is based 
on the adoption of higher-order Gauss-Hermite quadra- 
tures (linked to the discretization of the Boltzmann equa- 
tion), but its truncations yield at most asymptotic the- 
ories. Other approaches, which are based on 'ad hoc' 
modifications of the fluid equations (for example, intro- 
ducing compressibility and/or temperature effects [28j]). 
by definition cannot provide exact Navier-Stokes solvers. 

Another critical issue is related to the numerical sta- 
bility of LBM's 29], usually attributed to the violation of 
the condition of strict positivity (realizability condition) 
for the kinetic distribution function [2t| [3(|. Therefore, 
according to this viewpoint, a stability criterion should 
be achieved by imposing the existence of an H-theorem 
(for a review see 3l|). In an effort to improve the ef- 
ficiency of LBM numerical implementations and to cure 
these instabilities, there has been recently a renewed in- 
terest in the LB theory. Several approaches have been 
proposed. The first one involves the adoption of entropic 
LBM's (ELBM 0, [H [H, Q in which the equilibrium 



distribution satisfies also a maximum principle, defined 
with respect to a suitably defined entropy functional. 
However, usually these methods lead to non-polynomial 
equilibrium distribution functions which potentially re- 
sult in higher computational complexity [351 ] and less nu- 
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merical accuracy 36[ . Other approaches rely on the adop- 
tion of multiple relaxation times (3?], [HI . However the 
efficiency, of these methods is still in doubt. Therefore, 
the search for new [LB] models, overcoming these limita- 
tions, remains an important unsolved task. 



lb - Goals of the investigation 

The aim of this work is the development of an inverse 
kinetic theory for the incompressible Navier-Stokes equa- 
tions (INSE) which, besides realizing an exact Navier- 
Stokes (and Poisson) solver, overcomes some of the lim- 
itations of previous LBM's. Unlike Refs. [2, [3], where a 
continuous IKT was considered, here we construct a dis- 
crete theory based on the LB velocity-space discretiza- 
tion. In such a type of approach, the kinetic description 
is realized by a finite number of discrete distribution func- 
tions fi(r, t), for i = 0, k, each associated to a prescribed 
discrete constant velocity and defined everywhere in 
the existence domain of the fluid fields (the open set x / 
). The configuration space f2 is a bounded subset of the 
Euclidean space M 3 and the time interval / is a subset of 
R. The kinetic theory is obtained as in 0, 0| by introduc- 
ing an inverse kinetic equation (LB-IKE) which advances 
in time the distribution function and by properly defin- 
ing a correspondence principle, relating a set of velocity 
momenta with the relevant fluid fields. 

To achieve an IKT for INSE, however, also a proper 
treatment of the initial and boundary conditions, to be 
satisfied by the kinetic distribution function, must be in- 
cluded. In both cases, it is proven that they can be de- 
fined to be exactly consistent - at the same time - both 
with the hydrodynamic equations (which must hold also 
arbitrarily close to the boundary of the fluid domain) and 
with the prescription of the initial and Dirichlet bound- 
ary conditions set for the fluid fields. Remarkably, both 
the choice of the initial and equilibrium kinetic distri- 
bution functions and their functional class remain essen- 
tially arbitrary. In other words, provided suitable min- 
imal smoothness conditions are met by the kinetic dis- 
tributions function, for arbitrary initial and boundary 
kinetic distribution functions, the relevant moment equa- 
tions of the kinetic equation coincide identically with the 
relevant fluid equations. This includes the possibility 
of defining a LB-IKT in which the kinetic distribution 
function is not necessarily a Galilean invariant. 

This arbitrariness is reflected also in the choice of pos- 
sible "equilibrium" distribution functions, which remain 
essentially free in our theory, and can be made for exam- 
ple in order to achieve minimal algorithmic complexity. 
A possible solution corresponds to assume polynomial- 
type kinetic equilibria, as in the traditional asymptotic 
LBM's. These kinetic equilibria are well-known to be 
non- Galilean invariant with respect to arbitrary finite 
velocity translations. Nevertheless, as discussed in detail 



in Sec. 4, Subsection 4A, although the adoption of Galilei 
invariant kinetic distributions is in possible, this choice 
does not represent an obstacle for the formulation of a 
LB-IKT. Actually Galilean invariance need to be fulfilled 
only by the fluid equations. The same invariance prop- 
erty must be fulfilled only by the moment equations of 
the LB-IKT and not necessarily by the whole LB inverse 
kinetic equation (LB-IKE). 

Another significant development of the theory is the 
formal introduction of an cntropic principle, realized by 
a constant H-theorem, in order to assure the strict pos- 
itivity of the kinetic distribution function in the whole 
existence domain x I. The present entropic principle 
departs significantly from the literature. Unlike previ- 
ous entropic LBM's it is obtained without imposing any 
functional constraints on the class of the initial kinetic 
distribution functions. Namely without demanding the 
validity of a principle of entropy maximization (PEM, 
[39l ]) in a true functional sense on the form of the distri- 
bution function. Rather, it follows imposing a constraint 
only on a suitable set of extended fluid fields, in particu- 
lar the kinetic pressure pi(r,i).The latter is uniquely re- 
lated to the actual fluid pressure p(r, t) via the equation 
Pi(r,t) = p(r,t) + P {t), with P Q (t) > to be denoted 
as pseudo-pressure. The constant H-theorem is therefore 
obtained by suitably prescribing the function P (t) and 
implies the strict positivity. The same prescription as- 
sures that the entropy results maximal with respect in the 
class of the admissible kinetic pressures, i.e., it satisfies a 
principle of entropy maximization. Remarkably, since 
this property is not affected by the particular choice of 
the kinetic equilibrium, the H-thcorem applies also in the 
case of polynomial equilibria. We stress that the choice 
of the entropy functional remains essentially arbitrary, 
since no actual physical interpretation can be attached to 
it. For example, without loss of generality it can always 
be identified with the Gibbs-Shannon entropy. Even pre- 
scribing these additional properties, in principle infinite 
solutions exist to the problem. Hence, the freedom can 
be exploited to satisfy further requirements (for example, 
mathematical simplicity, minimal algorithmic complex- 
ity, etc.). Different possible realizations of the theory and 
comparisons with other CFD approaches are considered. 
The formulation of the inverse kinetic theory is also use- 
ful in order to determine the precise relationship between 
the LBM's and previous CFD schemes and in particular 
to obtain possible improved asymptotic LBM's with pre- 
scribed accuracy. As an application, we intend to con- 
struct asymptotic models which satisfy with prescribed 
accuracy the required fluid equations [INSE] and possi- 
bly extend also the range of validity of traditional LBM's. 
In particular, this permits to obtain asymptotic accuracy 
estimates of customary LB approaches. The scheme of 
presentation is as follows. In Sec. 2 the INSE problem 
is recalled and the definition of the extended fluid fields 
{V,pi} is presented. In Sec. 3 the basic assumptions 
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of previous asymptotic LBM's are recalled. In. Sec. 4 and 
5 the foundations of the new inverse kinetic theory are 
laid down and the integral LB inverse kinetic theory is 
presented, while in Sec. 6 the entropic theorem is proven 
to hold for the kinetic distribution function for properly 
defined kinetic pressure. Finally, in Sec. 7 various asymp- 
totic approximations are obtained for the inverse kinetic 
theory and comparisons are introduce with previous LB 
and CFD methods and in Sec. 8 the main conclusions 
are drawn. 



2 - THE INSE PROBLEM 



A prerequisite for the formulation of an inverse kinetic 
theory 0, @] providing a phase-space description of a clas- 
sical (or quantum) fluid is the proper identification of the 
complete set of fluid equations and of the related fluid 
fields. For a Newtonian incompressible fluid, referred to 
an arbitrary inertial reference frame, these are provided 
by the incompressible Navier-Stokes equations (INSE) for 
the fluid fields {p, V,p} 



(1) 
(2) 
(3) 

(4) 
(5) 

Equations ©-© are defined in a open connected set 
OCR 3 (defined as the subset of R 3 where p(r,t) > 0) 
with boundary SQ, while Eqs. ((4J and © apply on its 
closure il. Here the notation is standard. Thus, N is the 
NS operator 



V-V = 0, 
NV = 0, 

p{v,t) = p . 

There are supplemented by the inequalities 

p(r,t) > 0, 
Po > 0. 



f- M V 2 V, 



(6) 



with^ 



+ V • V the convective derivative, f denotes 
a suitably smooth volume force density acting on the fluid 
element and p = vp Q > is the constant fluid viscosity. 
In particular we shall assume that f can be represented 
in the form 



-V$( 



■fi(r,t) 



where we have separated the conservative V$(r) and the 
non-conservative f\ parts of the force. Equations (HI)-© 
are assumed to admit a strong solution in Q x J, with 
I Cla possibly bounded time interval. By assumption 
{p, V,p} are continuous in the closure O. Hence if in £1 x I, 
f is at least C^ 1,0 \flxl), it follows necessarily that {V,p} 
must be at least C^ 2,1 \Vl x 7). In the sequel we shall 
impose on {V,p} the initial conditions 



V(r,t„) - V Q (r), 
p(r,t a ) = p {r). 



(7) 



Furthermore, for greater mathematical simplicity, here 
we shall impose Dirichlet boundary conditions on 5tt 



V(-,t)| tfn = V w (;t)\ sn 
K->*)l«i = Pw(;t)\ sn . 



(8) 



Eqs.© and ©-© define the initial-boundary value 
problem associated to the reduced INSE (reduced INSE 
problem). It is important to stress that the previous 
problem can also formulated in an equivalent way by re- 
placing the fluid pressure p(r, t) with a function pi (r, t) 
(denoted kinetic pressure) of the form 



pi(r,t) = P +p(r,t), 



(9) 



where P — P Q (t) is prescribed (but arbitrary) real func- 
tion of time and is at least P a (t) e C^(I). {V,pi} will 
be denoted hereon as extended fluid fields and P (t) will 
be denoted as pseudo-pressure. 



3 - ASYMPTOTIC LBM'S 
3A - Basic assumptions 

As is well known, all LB methods are based on a dis- 
crete kinetic theory, using a so-called lattice Boltzmann 
velocity discretization of phase-space (LB discretization) . 
This involves the definition of a kinetic distribution func- 
tion /, which can only take the values belonging to a 
finite discrete set {/j(r, t),i — 0,k} (discrete kinetic dis- 
tribution functions) . In particular, it is assumed that the 
functions fi, for % = 0, k, are associated to a discrete set 
of k + 1 different "velocities" {a^, i — 0, k) . Each is an 
'a priori' prescribed constant vector spanning the vector 
space W 1 (with n — 2 or 3 respectively for the treatment 
of two- and three-dimensional fluid dynamics), and each 
fi(r,t) is represented by a suitably smooth real function 
which is defined and continuous in f2 x I and in particular 
is at least C^ k ' j \Q, x I) with k > 3. 

The crucial aspect which characterizes customary LB 



ine crucial aspect wmcn cnaracterizes customary lb 
approaches 0, S S3, E3, E3, E3, E3, S3, E^l involves the 
construction of kinetic models which allow a finite sound 
speed in the fluid and hence are based on the assumption 
of a (weak) compressibility of the same fluid. This is 
realized by assuming that the evolution equation (kinetic 
equation) for the discrete distributions fi(r,t) (i — l,k), 
depends at least one ( or more ) infinitesimal ( asymptotic ) 
parameters (see below). Such approaches are therefore 
denoted as asymptotic LBM's. They are characterized 
by a suitable set of assumptions, which typically include: 

1. LB assumption #1: discrete kinetic equation and 
correspondence principle: the first assumption con- 
cerns the definition of an appropriate evolution 
equation for each /i(r, t) which must hold (together 
with all its moment equations) in the whole open 
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set f2 x /. In customary LB approaches it takes the 
form of the so-called LB-BGK equation HHUm 



L(i)fi = Oi(/i), 



(10) 



where i = 0, k. Here L(j\ is a suitable streaming 
operator, 



Wi) = -vcVi - n q ) 



(ii) 



(with v c > a constant collision frequency) is 
known as BKG collision operator (after Bhatba- 
gar, Gross and Krook 143] ) and fl q is an "equi- 
librium" distribution to be suitably defined. In 
customary LBM's it is implicitly assumed that 
the solution of Eq. (flO|) . subject to suitable initial 
and boundary conditions exists and is unique in 
the functional class indicated above. In partic- 
ular, usually L(j\ is either identified with the fi- 
nite difference streaming operator (see for example 
[I El, H, El), i.e., L {i) fi(r,t) = L FD(i) /i(r,t) = 
^ [fi(r + a,Ai, t + At) ~ /<(r, t)] or with the dif- 
ferential streaming operator (see for instance [171 



L(i) = L D{l) = — 



d_ 
dr 



(12) 



Here the notation is standard. In particular, in the 
case of the operator Lpntiji <-^t and cAt = L Q are 
appropriate parameters which define respectively 
the characteristic time- and length- scales associ- 
ated to the LBM time and spatial discretizations. 
A common element to all LBM's is the assump- 
tion that all relevant fluid fields can be identified, 
at least in some approximate sense, with appro- 
priate momenta of the discrete kinetic distribu- 
tion function {correspondence principle). In par- 
ticular, for neutral and isothermal incompressible 
fluids, for which the fluid fields are provided re- 
spectively by the velocity and pressure fluid fields 
{Y j (T,t),j = l,4} = {V(r,t),p(r,t)}, it is as- 
sumed that they are identified with a suitable set 
of discrete velocity momenta (for j — 1,4) 



Yj(r,t) = X Tl (r,t)Mr,t), 



(13) 



i=0, k 



where Xji(r, t) (with i = 0, k and j = l,k) are ap- 
propriate, smooth real weight functions. In the 
literature several examples of correspondence prin- 
ciples are provided, a particular case b eing provided 



by the so-called D2Q9 (V,p)-scheme [4 



p(r,t)= C 2 ^ / ^c 2 ^ / ^, 



(14) 



i=0, k 



-0,k 



V(r,t) = i-E *** = - E a ^ eq) > ( 15 ) 



where k — 8 and c = min{|a i | > 0,i = 0, k} is a 
characteristic parameter of the kinetic model to be 
interpreted as test particle velocity. In customary 
LBM's the parameter c s = (with D the dimen- 
sion of the set Q) is interpreted as sound speed of 
the fluid. In order that the momenta (fLT| and (fT5|) 
recover (in some suitable approximate sense) INSE 
, however, appropriate subsidiary conditions must 
be met. 

2. LB assumption #2: Constraints and asymptotic 
conditions: these are based on the introduc- 
tion of a dimensionless parameter e, to be consid- 
ered infinitesimal, in terms of which all relevant 
parameters can be ordered. In particular, it is 
required that the following asymptotic orderings 
fvA liol l4l| apply respectively to the fluid fields 
po, V(r, t),p(r, t), the kinematic viscosity v — p/p 
and Reynolds number R e = LV/u: 

p o ,V(v,t),p{r,t)~o(e ), (16) 

(17) 



v = — [l + o{e)]~o{s a »), 



l/o{e a "), 



(18) 



where aji > 0. Here we stress that the position 
for v holds in the case of D2Q9 only, while the 
generalization to 3D and other LB discretizations, 
is straightforward. Furthermore, the velocity c and 
collision frequency v c are ordered so that 



c~l/o(^), 
Vc ~l/o{e a ») 
-^~o(e a ), 



(19) 
(20) 
(21) 



with a = a v — a c > 0; the characteristic length and 
time scales, L Q = cAt and At for the spatial and 
time discretization are assumed to scale as 

° At L ° - (22) 



L L 



O £ 



(23) 



with a>t,C(L > 0. Here L and T are the (smallest) 
characteristic length and time scales, respectively 
for spatial and time variations of V(r, t) and p(r.i). 
Imposing also that results infinitesimal at least 
of order 

7^-0(6") 
1 V c 

it follows that it must be also at — oil > 0. These 
assumptions imply necessarily that the dimension- 
less parameter M e ^ = — (Mach number) must be 
ordered as 



Po 



1=1, fc 



Po 



1=1, fc 



M eff ~ 0{e a ") 
{small Mach-number expansion) 



(24) 
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3. LB assumption #3: Chapman- Enskog expansion - 
Kinetic initial conditions, relaxation conditions: it 
is assumed that the kinetic distribution function 
fi(r,t) admits a convergent Chapman-Enskog ex- 
pansion of the form 



r2 ;(2) 



(25) 



where S = e a and the functions (j £ N) 
are assumed smooth functions of the form (multi- 
scale expansion) fi(r a , ri, ra, ..t , t\, ti, ..), where 
r n = 5 n r, t n = 5 n t and n G N. In typical LBM's 
the parameter S is usually identified with e (which 
requires letting a = 1), while the Chapman-Enskog 
expansion is usually required to hold at least up to 
order o(S 2 ). In addition the initial conditions 



fi(r,t ) = f° q (r,t ) 



(26) 



(for i = 0, k) are imposed in the closure of the fluid 
domain f2. It is well known [46[ that this position 
generally (i.e., for non-stationary fluid fields), im- 
plies the violation of the Chapman-Enskog expan- 
sion close to t = t a , since the approximate fluid 
equations are recovered only letting Sf^ +S 2 f^ ^ 
0, i.e., assuming that the kinetic distribution func- 
tion has relaxed to the Chapman-Enskog form ([25]) . 
This implies a numerical error (in the evaluation of 
the correct fluid fields) which can be overcome only 
discarding the first few time steps in the numerical 
simulation. 

4. LB assumption #5: Equilibrium kinetic distribu- 
tion: a possible realization for the equilibrium dis- 
tributions f^ q (i = 0, k) is given by a polynomial of 
second degree in the fluid velocity [4J| 



(27) 



+W % p 



a, V 



a, V 



II! 

2~?~ 



Here, without loss of generality, the case of the 
D2Q9 LB discretization will be considered, with mj 
and ELi (for i = 0, 8) denoting prescribed dimension- 
less constant weights and discrete velocities. Notice 
that, by definition, f^ q is not a Galilei scalar. Nev- 
ertheless, it can be considered approximately in- 
variant, at least with respect to low- velocity trans- 
lations which do not violate the low-Mach number 
assumption (j24|) . 

5. LB assumption #6: Kinetic boundary conditions: 
They are specified by suitably prescribing the form 
of the inco ming distribution function at the bound- 
ary 5fi . El m, 0, M, 0, 001IF 



However, this position is not generally 



consistent with the Chapman-Enskog solution (|25|) 
(see related discussion in Appendix A). As a con- 
sequence violations of the hydrodynamic equations 
may be expected sufficiently close to the boundary, 
a fact which may be only alleviated (but not com- 
pletely eliminated) by adopting suitable grid refine- 
ments near the boundary. An additional potential 
difficulty is related to the condition of strict posi- 
tivity of the kinetic distribution function [57J which 
is not easily incorporated into the no-slip boundary 



conditions 



3B - Computational complexity of asymptotic 
LBM's 



The requirements posed by the validity of these hy- 
potheses may strongly influence the computational com- 
plexity of asymptotic LBM's which is usually associated 
to the total number of "logical" operations which must 
be performed during a prescribed time interval. There- 
fore, a critical parameter of numerical simulation meth- 
ods is their discretization time scale At. This is - in turn 
- related to the Courant number Nq = ^7^, where V 
and L . denote respectively the sup of the magnitude of 
the fluid velocity and the amplitudes of the spatial dis- 
cretization. As is well known "optimal" CFD simulation 
methods typically allow L a ~ L and a definition of the 
time step At = At 0p t such that N c ~ VA ^°" t ~ 1. In- 
stead, for usual LBM's satisfying the low-M 6 -'-' assump- 
tion (|24[) . the Courant number is very small since it re- 
sults N c = Af e// % ~ 0{e a )^. Tnis means that their 
discretization time scale of At is much smaller than Ato P t 
and reads 



At 



M eff — Ai 
L 



Opt- 



(28) 



In addition, depending on the accuracy of the numeri- 
cal algorithms adopted for the construction of the dis- 
crete kinetic distribution function, also the ratio re- 
sults infinitesimal in the sense ~ o(e aL ), with suitable 
oil > 0. Finally, we stress that LB approaches based 
on the adoption of the finite-difference streaming opera- 
tor Lfd^) are usually only accurate to order o(At 2 ). For 
them, therefore, the requirement placed by Eq. (|2"5|) might 
be even stronger. This implies that traditional LBM's 
may involve a vastly larger computation time than that 
afforded by more efficient numerical methods. 



NEW LB INVERSE KINETIC THEORY 
(LB-IKT) 



42] 



con- 



A basic issue in LB approaches [8J, |LJJ, [li 
cerns the choice of the functional class of the discrete 
kinetic distribution functions /$ (i — 0, k) as well as the 
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related definition of the equilibrium discrete distribution 
function /■ [which appears in the BGK collision opera- 
tor; see Eq. fTTj) ]. This refers in particular to their trans- 
formation properties with respect to arbitrary Galilean 
transformations, and specifically to their Galilei invari- 
ance with respect to velocity translations with constant 
velocity. 

In statistical mechanics it is well known that the ki- 
netic distribution function is usually assumed to be a 
Galilean scalar. The same assumption can, in principle, 
be adopted also for LB models. However, the kinetic 
distribution functions fi and f^ q do not necessarily re- 
quire a physical interpretation of this type. In the se- 
quel we show that for a discrete inverse kinetic theory it 
is sufficient that fi and f^ q be so defined that the mo- 
ment equations coincide with the fluid equations (which 
by definition are Galilei covariant). It is sufficient to de- 
mand that both fi and f^ q are identified with a ordinary 
scalars with respect to the group of rotation in IR 2 , while 
they need not be necessarily invariant with respect to 
arbitrary velocity translations. This means that fi is in- 
variant only for a particular subset of inertial reference 
frames. For example for a fluid which at the initial time 
moves locally with constant velocity an element of this 
set can be identified with the inertial frame which in the 
same position is locally co-moving with the fluid. 

The adoption of non-translationally invariant discrete 
distributions fi is actually already well known in LBM 
and results convenient for its simplicity. This means, 
manifestly, that in general no obvious physical interpre- 
tation can be attached to the other momenta of the dis- 
crete kinetic distribution function. As a consequence, 
the very definition of the concept of statistical entropy to 
be associated to the f[s is essentially arbitrary, as well as 
the related principle of entropy maximization, typically 
used for the determination of the equilibrium distribution 
function f^ q . Several authors, nevertheless, have investi- 
gated the adoption of possible alternative formulations, 
which are based on suitable definitions of the entropy 
functional and/or the requirement of approximate or ex- 
act Galilei invariance (see for example [29|, 32, 62l|). 



should yield exactly the complete set of fluid equations 
and which, contrary to customary kinetic approaches in 
CFD (in particular LB methods), should not depend on 
asymptotic parameters. This implies that the inverse ki- 
netic theory must also satisfy an exact closure condition. 
As a further condition, we require that the fluid equa- 
tions are fulfilled independently of the initial conditions 
for the kinetic distribution function (to be properly set) 
and should hold for arbitrary fluid fields. The latter re- 
quirement is necessary since we must expect that the 
validity of the inverse kinetic theory should not be lim- 
ited to a subset of possible fluid motions nor depend on 
special assumptions, like a prescribed range of Reynolds 
numbers. In principle a phase-space theory, yielding an 
inverse kinetic theory, may be conveniently set in terms of 
a quasi-probability, denoted as kinetic distribution func- 
tion, /(x,i). A particular case of interest (investigated in 
Refs.[2j, y]) refers to the case in which /(x, i) can actu- 
ally be identified with a phase-space probability density. 
In the sequel we address both cases, showing that, to a 
certain extent, in both cases the formulation of a generic 
IKT can actually be treated in a similar fashion. This 
requires the introduction of an appropriate set of consti- 
tutive assumptions (or axioms). These concern in par- 
ticular the definitions of the kinetic equation - denoted 
as inverse kinetic equation (IKE) - which advances in 
time /(x, t) and of the velocity momenta to be identified 
with the relevant fluid fields {correspondence principle). 
However, further assumptions, such as those involving 
the regularity conditions for /(x, t) and the prescription 
of its initial and boundary conditions must clearly be 
added. The concept [of IKT] can be easily extended to 
the case in which the kinetic distribution function takes 
on only discrete values in velocity space. In the sequel 
we consider for definiteness the case of the so-called LB 
discretization, whereby - for each (r,t) £ Q x I - the 
kinetic distribution function is discrete, and in particu- 
lar admits a finite set of discrete values fi(r,t) € K, for 
i = 0, k, each one corresponding to a prescribed constant 
discrete velocity 6 K 3 for i = 0, k. 



4A - Foundations of LB-IKT 

As previously indicated, there are several important 
motivations for seeking an exact solver based on LBM. 
The lack of a theory of this type represents in fact a 
weak point of LB theory. Besides being a still unsolved 
theoretical issue, the problem is relevant in order to de- 
termine the exact relationship between the LBM's and 
traditional CFD schemes based on the direct discretiza- 
tion of the Navier-Stokes equations. Following ideas re- 
cently developed @,[3,I3,[EQ> we show that such a theory 
can be formulated by means of an inverse kinetic theory 
(IKT) with discrete velocities. By definition such an IKT 



4B - Constitutive assumptions 

Let us now introduce the constitutive assumptions (ax- 
ioms) set for the construction of a LB-IKT for INSE, 
whose form is suggested by the analogous continuous 
inverse kinetic theory 0, Q. The axioms, define the 
"generic" form of the discrete kinetic equation, its func- 
tional setting, the momenta of the kinetic distribution 
function and their initial and boundary conditions, are 
the following ones: 
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Axiom I - LB-IKE and functional setting. 



initial conditions 



Let us require that the extended fluid fields {V,Pi} 
are strong solutions of INSE, with initial and boundary 
conditions (JZ}-® and that the pseudo pressure p (t) is 
an arbitrary, suitably smooth, real function. In particu- 
lar we impose that the fluid fields and the volume force 
belong to the minimal functional setting: 

Pi,$eC (2a) (ft x /), 
VeC (3 ' 1} (fi x /), (29) 
fieC (1 ' 0) (^ x I). 
We assume that in the set O x / the following equation 



fi(r,t ) = foi(r,t ) 



(33) 



(30) 



[LB inverse kinetic equation (LB-IKE)] is satisfied iden- 
tically by the discrete kinetic distributions fi(r,t) for 
i = 0, k. Here Oi(fi) and L D ^ are respectively the BGK 
and the differential streaming and operators [Eqs. (fTT|) 
and (fl~2|) ]. while Si is a source term to be defined. We 
require that KB- IKE is defined in the set f2 x J, so that 
Oi(fi) and Si are at least that C^(0 x /) and contin- 
uous in fi x J. Moreover Oi(fi), defined by Eq. (fTTj) . is 
considered for generality and will be useful for compar- 
isons with customary LB approaches. We remark that 
the choice of the equilibrium kinetic distribution ff q in 
the BGK operator remains completely arbitrary. We 
assume furthermore that in terms of fi the fluid fields 
{V,pi} are determined by means of functionals of the 
form Mx [fi] — ^2 Xjfi (denoted as discrete velocity 

' i=0,8 

momenta). For X = X\,X 2 (with X\ — c 2 ,X 2 — 
these are related to the fluid fields by means of the equa- 
tions (correspondence principle) 



Pi 



(r,i)-$(r)=c 2 ^/. i =c 2 ^/r, (31) 



i=0,8 



V(r,t)=l E aifi = iE > 



Po 



Po 



(32) 



where c = min{|ai| , i = 1,8} is the test particle veloc- 
ity and fl q is defined by Eq. (j2"Tj) but with the kinetic 
pressure p\ that replaces the fluid pressure p adopted 
previously 44 1. These equations are assumed to hold 



identically in the set fi x I and by assumption, fi and 
fl q belong to the same functional class of real functions 
defined so that the extended fluid fields belong to the 
minimal functional setting (|29p. Moreover, without loss 
of generality, we consider the D2Q9 LB discretization. 



Axiom II - Kinetic initial and boundary conditions. 

The discrete kinetic distribution function satisfies, for 
% = 0, k and for all r belonging to the closure fi, the 



where / i(r,i ) (for i = 0, k) is a initial distribution func- 
tion defined in such a way to satisfy in the same set the 
initial conditions for the fluid fields 

p Xo (r) = P {t ) + Po (r) - 4>(r) = (34) 
= c 2 E / oi (r), 



V (r) = - V aJ OT (r) 

Po ^ 



(35) 



i=l,8 



To define the analogous kinetic boundary conditions on 
5fi, let us assume that <5fi is a smooth, possibly moving, 
surface. Let us introduce the velocity of the point of the 
boundary determined by the position vector r w £ <5fi, de- 
fined by V w (r w (t),t) = 4jX w {t) and denote by n(r tu ,t) 
the outward normal unit vector, orthogonal to the bound- 
ary 80 at the point r w . Let us denote by f^(r w ,t) 
and f\ ' (r w , t) the kinetic distributions which carry the 
discrete velocities for which there results respectively 
(a* — V w ) ■ n(r w ,t) > (outgoing- velocity distributions) 
and (an — ~V W ) ■ n(r w ,t) < (incoming- velocity distribu- 
tions) and which are identically zero otherwise. We as- 
sume for definiteness that both sets, for which |aj| > 0, 
are non empty (which requires that the parameter c 
be suitably defined so that c > |V m |). The bound- 
ary conditions are obtained by prescribing the incom- 
ing kinetic distribution fj (r w ,t), i.e., imposing (for all 
(r w ,t) £ SO x I) 



ft\T w ,t) = f { j\r w ,t). 



(36) 



Here f^~\v Wl t) are suitable functions, to be assumed 
non- vanishing and defined only for incoming discrete ve- 
locities for which (a, — V w )-n(r w ,t) < 0. Manifestly, the 
functions \r w ,t) (i = 0, k) must be defined so that 
the Dirichlet boundary conditions for the fluid fields are 
identically fulfilled, namely there results 



Piw(r w ,t) = P (t) +p w (r w ,t) - $(r) 
E {f i o ;\r w ,t) + ft ) (r w ,t)} ) 



i=0,k 



V w (r w ,t) = 

= ~ E ^{f { o ;\r w ,t) + ri +) (r w ,t)\. 



(37) 



(38) 



Here, again, the functions / i(r) and fof\^w,t) (for i — 
0, k) must be assumed suitably smooth. A particular case 
is obtained imposing identically for i = 0, k 



foi(*,to) 

f { ot\r w ,t) 



ifMo), 
f: 9 (r w ,t), 



(39) 
(40) 
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where the identification with f„f\r w ,t) and f^i'{r w ,t) 
is intended respectively in the subsets aj-n(r„,, t) > and 
&i -n(r„,, t) < 0. Finally, we notice that in case Neumann 
boundary conditions are imposed on the fluid pressure, 
Eq. ([57|) still holds provided p w (r w ,t) is intended as a 
calculated value. 



f(-), 



Axiom III - Moment equations. 

If fi(r,t), for i = 0, k, are arbitrary solutions of LB- 
IKE [Eg. ([30]) ] which satisfy Axioms I and II validity of 
Axioms I and II, we assume that the moment equations of 
the same LB-IKE, evaluated in terms of the moment op- 
erators Mx [•] = Xj-, with j = 1,2, coincide iden- 

3 i=0,8 

tically with INSE, namely that there results identically 
[for all (r, t) € Q X I] 

M Xl [L i fi-n i (f i )-S i ]=V-V = 0, (41) 



M X2 \UU - fli(fi) - S,} = NV = 0. (42) 



Axiom IV - Source term. 

The source term is required to depend on a finite num- 
ber of momenta of the distribution function. It is as- 
sumed that these include, at most, the extended fluid 
fields {V,pi} and the kinetic tensor pressure 



H = 3j2f iaiai -p VV. 



(43) 



• Furthermore, we also normally require (except 
for the LB-IKT described in Appendix B) that 
Si(r,t) results independent of /? 9 (r,i), f i(r) and 
fwi(r w ,t) (for i = Q,k). 

Although, the implications will made clear in the fol- 
lowing sections, it is manifest that these axioms do not 
specify uniquely the form (and functional class) of the 
equilibrium kinetic distribution function / 4 ei? (r,i), nor 
of the initial and boundary kinetic distribution func- 
tions (J33J) , d36j) . Thus, both f- q (r,t), / oi (r,t ) and the 
related distribution they still remain in principle com- 
pletely arbitrary. Nevertheless, by construction, the 
initial and (Dirichlet) boundary conditions for the fluid 
fields are satisfied identically. In the sequel we show that 
these axioms define a (non-empty) family of parameter- 
dependent LB-IKT 's, depending on two constant free pa- 
rameters v c , c > and one arbitrary real function P (t). 
The examples considered are reported respectively in the 
following Sec. 5,6 and in the Appendix B. 



5 - A POSSIBLE REALIZATION: THE 
INTEGRAL LB-IKT 

We now show that, for arbitrary choices of the distri- 
butions fi(v,t) and fi 9 (r,t) which fulfill axioms I-IV, an 
explicit (and non-unique) realization of the LB-IKT can 
actually be obtained. We prove, in particular, that a pos- 
sible realization of the discrete inverse kinetic theory, to 
be denoted as integral LB-IKT, is provided by the source 
term 



Wi 
72" 



^ - a< • (fi-^v 2 v - v • n + Vp) 



(44) 

= Si, 



where is denoted as first pressure term. Holds, in 

fact, the following theorem. 



Theorem 1 - Integral LB-IKT 

In validity of axioms I-IV the following statements 
hold. For an arbitrary particular solution fi and for ar- 
bitrary extended fluid fields: 

A) if fi is a solution of LB-IKE [Eg. i30\) l the moment 
equations coincide identically with INSE in the set Six I; 

B) the initial conditions and the (Dirichlet) boundary 
conditions for the fluid fields are satisfied identically; 

C) in validity of axiom IV the source term Si is non- 
uniquely defined by Eg. {44\ l- 

Proof 

A ) We notice that by definition there results identically 



i=0 



Si 



1 dpi 

c 2 dt 



(45) 



(46) 



i=0 



= -- [f-/iV 2 v-v-n+ vp\ 

On the other hand, by construction (Axiom I) fi (i — 
1, k) is defined so that there results identically Y^,i=o ^ = 
and 53i=o a »^i = 0. Hence the momenta M Xl , Mx 2 of 
LB-IKE deliver respectively 



V • a*/* = 



(47) 



3 di a ^ + P° Y ■ VV + V ^ + f -MV 2 V = (48) 



»=i,e 
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where the fluid fields V,pi are defined by Eqs. (f3T]) . (l32p . 
Hence Eqs. (|47|) and (|48|) coincide respectively with the 
isochoricity and Navier-Stokes equations [((T|) and @]. 
As a consequence, is a particular solution of LB-IKE 
iff the fluid fields {V,pi} are strong solutions of INSE. 

B) Initial and boundary conditions for the fluid fields 
are satisfied identically by construction thanks to Axiom 
II. 

C) However, even prescribing v c , c > and the real 
function P (t), the functional form of the equation can- 
not be unique The non uniqueness of the functional form 
of the source term Si(r,t) is assumed to be indepen- 
dent of fi 9 (r,t) [and hence of Eq. ((30|) ] is obvious. In 
fact, let us assume that Si is a particular solution for 
the source term which satisfies the previous axioms I- 
IV. Then, it is always possible to add to Si arbitrary 
terms of the form Si + 5 Si, with 5 Si ^ which depends 
only on the momenta indicated above, and gives van- 
ishing contributions to the first two moment equations, 
namely M Xj [SSi] = £ XjSSi = 0, with j = 1,2. To 

i=0,8 

prove the non-uniqueness of the source term Si, it is suf- 
ficient to notice that, for example, any term of the form 
SSi = f — 1^ F(r,t), with F(r,t) an arbitrary real 

function (to be assumed, thanks to Axiom IV, a linear 
function of the fluid velocity), gives vanishing contribu- 
tions to the momenta Mx t , Mx 2 ■ Hence Si is non-unique. 

The implications of the theorem are straightforward. 
First, manifestly, it holds also in the case in which the 
BGK operator vanishes identically. This occurs letting 
v c = in the whole domain Q x I. Hence the inverse 
kinetic equation holds independently of the specific defi- 
nition of /? 9 (r,i). 

An interesting feature of the present approach lies in 
the choice of the boundary condition adopted for /j(r,i), 
which is different from that usually adopted in LBM's 
[see for example [14j for a review on the subject]. In par- 
ticular, the choice adopted is the simplest permitting to 
fulfill the Dirichlet boundary conditions [imposed on the 
fluid fields]. This is obtained prescribing the functional 
form of fi(r,t) on the boundary of the fluid domain (SQ), 
which is identified with a function f i(r,t). 

Second, the functional class of fi(r,t), f^ q (r,t) and of 
f i(r,t) remains essentially arbitrary. Thus, in particu- 
lar, the initial and boundary conditions, specified by the 
same function / j(r, t), can be defined imposing the po- 
sitions ([39]) . (|40 | . As further basic consequence, fP{r,t) 
and fi(r,t) need not necessarily be Galilei-invariant (in 
particular they may not be invariant with respect to ve- 
locity translations), although the fluid equations must 
be necessarily fully Galilei-covariant. As a consequence 
it is always possible to select /? 9 (r,t) and f i(r,t) based 
on convenience and mathematical simplicity. Thus, be- 
sides distributions which are Galilei invariant and sat- 
isfy a principle of maximum entropy (see for example 



[2 21 . l30l . |32j, [3J, [60|, |6l|), it is always possible to iden- 
tify them [i.e., f^ q (r,t), f i(r,t)] with a non-Galilean in- 
variant polynomial distribution of the type (|27j) [mani- 
festly, to be exactly Galilei-invariant each f^ q (r,t) should 
depend on velocity only via the relative velocity = 
a*-V]. 

We mention that the non-uniqueness of the source term 
Si can be exploited also by imposing that f? q (r,t) re- 
sults a particular solution of the inverse kinetic equation 
Eg. ([30)1 and there results also / i(r, t) = f^ q (r,t). In Ap- 
pendix B we report the extension of THM.l which is ob- 
tained by identifying again f^ q (v,t) with the polynomial 
distribution (|27|) . 



THE ENTROPIC PRINCIPLE - CONDITION 
OF POSITIVITY OF THE KINETIC 
DISTRIBUTION FUNCTION 



A fundamental limitation of the standard LB ap- 
proaches is their difficulty to attain low viscosities, due to 
the appearance of numerical instabilities [3] ■ In numeri- 
cal simulations based on customary LB approaches large 
Reynolds numbers is usually achieved by increasing nu- 
merical accuracy, in particular strongly reducing the time 
step and the grid size of the spatial discretization (both 
of which can be realized by means of numerical schemes 
with adaptive time-step and using grid refinements). 
Hence, the control [and possible inhibition] of numerical 
instabilities is achieved at the expense of computational 
efficiency. This obstacle is only partially alleviated by 
approaches based on ELBM [llllESSSIill. Such 
methods are based on the hypothesis of fulfilling an H- 
theorem, i.e., of satisfying in the whole domain £1 x / 
the condition of strict positivity for the discrete kinetic 
distribution functions. This requirement is considered, 
by several authors (see for example 2(| 2!^ 6^1), an es- 
sential prerequisite to achieve numerical stability in LB 
simulations. However, the numerical implementation of 
ELBM typically induce a substantial complication of the 
original algorithm, or require a cumbersome fine-tuning 
of adjustable parameters 22l. 37|. 



6A - The constant entropy principle and PEM 

A basic aspect of the IKT's here developed is the possi- 
bility of fulfilling identically the strict positivity require- 
ment by means of a suitable H-theorem which provides 
also a maximum entropy principle. In particular, in this 
Section, extending the results of THM.l and 2, we intend 
to prove that a constant H-theorem can be established 
both for the integral and differential LB- IKT's defined 
above. The H-theorem can be reached by imposing for 
the Gibbs-Shannon entropy functional the requirement 
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that for alH e / there results 



d_ 

at 



s(f) 



d_ 

"dt 



fd 3 r ^ /iln(/iM)=0, (49) 



i=0,8 



which implies that is necessarily maximal in a suit- 
able functional set {/} . The result can be stated as fol- 
lows: 



Theorem 2 - Constant H-theorem 

In validity of THM.l, let us assume that: 

1 ) the configuration domain f2 is bounded; 

2) at time t the discrete kinetic distribution functions 
fi, for i = 0, 8, are all strictly positive in the set f2. 

Then the following statements hold: 

A) by suitable definition of the pseudo pressure 
Po{t)-, the Gibbs-Shannon entropy functional S(f) = 
— J d 3 r ^2 fi^{fi/ w i) can be set to be constant in the 

n i=0,8 

whole time interval I. This holds provided the pseudo- 
pressure P (t) satisfies the differential equation 

^/d 3 r£5(l+log/0= (50) 
f 8 

= / d 3 rJ2(^-^fi-Si) (1 + log/i), 



where % = # + f-ff ; 

B) if the entropy functional S(f) = 
— J d 3 r h^ n {fil w i) * s constant in the whole 

Q i=0,8 

time interval I the discrete kinetic distribution functions 
fi are all strictly positive in the whole set Q x I; 

C) an arbitrary solution of LB-IKE [Ea.f30\)J which 
satisfies the requirement A) is extremal in a suitable func- 
tional class and maximizes the Gibbs-Shannon entropy . 

Proof: 

A) Invoking Eg. ([50)1 . there results 



dS(t) 
dt 



o i=0 



r 8 

= / d 3 rJ2(*i- V fi ~ Si) (1 + log 

J a n 



i=0 



where Si is the source term, provided by Eq. (l44|) . By 
direct substitution it follows the thesis. 

B) If Eg. ([50)) holds identically in there results W G 
J, S (t) = S (io) j which implies the strict positivity of /j, 
for all i = 0, 8. 

C) Let us introduce the functional class 

{/ + aSf} = {fi = fi(t) + aSfi(t),i = 0, 8} , (52) 



where a is a finite real parameter and the syn- 
chronous variation Sfi(t) is defined Sfi(t) — dfiit) = 
dt. Introducing the synchronous variation of the en- 

with ijj(a) = 



dt 



tropy, defined by SS (t) 
S(f + aSf) , it follows 



a=0 



SS (t) = dt 



dS{t) 
dt ■ 



(53) 



Since in validity of Eg. ([50)1 there results 



dS(t) _ 

at 



which in view of Eg. ([53)) implies also SS (t) = 0. It is im- 
mediately follows that there results necessarily S 2 S (t) < 
0, i.e., S (t) is maximal. Therefore, the kinetic distribu- 
tion function which satisfies IKE fEq.([30))] is extremal in 
the functional class of variations (|52p and maximizes the 
Gibbs-Shannon entropy functional. 



6B - Implications 

In view of statement B, THM.2 warrants the strict pos- 
itivity of the discrete distribution functions /, (i — 0, 8) 
only in the open set x J, while nothing can be said 
regarding their behavior on the boundary SQ (on which 
fi might locally vanish). However, since the inverse ki- 
netic equation actually holds only in the open set £1x1, 
this does not affect the validity of the result. While the 
precise cause of the numerical instability of LBM's is still 
unknown, the strict positivity of the distribution function 
is usually considered important for the stability of the nu- 
merical solution [29|, [33| ■ It must be stressed that the nu- 
merical implementation of the condition of constant en- 
tropy Eq. ()50p should be straightforward, without involv- 
ing a significant computational overhead for LB simula- 
tions. Therefore it might represent a convenient scheme 
to be adopted also for customary LB methods. 



7 - ASYMPTOTIC APPROXIMATIONS AND 
COMPARISONS WITH PREVIOUS CFD 
METHODS 

A basic issue is the relationship with previous CFD nu- 
merical methods, particularly asymptotic LBM's. Here 
we consider, for definiteness, only the case of the inte- 
gral LB-IKT introduced in Sec. 5. Another motivation is 
the possibility of constructing new improved asymptotic 
models, which satisfy with prescribed accuracy the re- 
quired fluid equations [INSE] , of extending the range of 
validity of traditional LBM's and fulfilling also the en- 
tropic principle (see Sec. 6). The analysis is useful in 
particular to establish on rigorous grounds the consis- 
tency of previous LBM's. The connection [with previ- 
ous LBM's] can be reached by introducing appropriate 
asymptotic approximations for the IKT's, obtained by 
assuming that suitable parameters which characterize the 
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IKT's are infinitesimal (or infinite) (asymptotic parame- 
ters). A further interesting feature is the possibility of 
constructing in principle a class of new asymptotic LBM's 
with prescribed accuracy , i.e., in which the distribution 
function (and the corresponding momenta) can be de- 
termined with predetermined accuracy in terms of per- 
turbative expansions in the relevant asymptotic parame- 
ters. Besides recovering the traditional low-Mach number 
LBM's fTJ ' 



parameters 



21 



401 ] . which satisfy the isochoricity condi- 
tion only in an asymptotic sense and are closely related to 
the Chorin artificial compressibility method, it is possible 
to obtain an improved asymptotic LBM's which satisfy 
exactly the same equation. 

We first notice that the present IKT is characterized 
by the arbitrary positive parameters v Cl c and the initial 
value P (t ), which enter respectively in the definition 
of the BGK operator [see Ijlip], the velocity momenta 
and equilibrium distribution function f^ q . Both c and 
P {t ) must be assumed strictly positive, while, to assure 
the validity of THM.2, P (t ) must be defined so that 
(for all i = 0,8) f- q (r,t ) > in the closure Q. Thanks 
to THM.l.and 2 the new theory is manifestly valid for 
arbitrary finite value of these parameters. This means 
that they hold also assuming 



(54) 

(55) 
(56) 



o(e a »)' 
1 



o(e<*<) ' 
Po(t ) ~ o(e°), 



where e denotes a strictly positive real infinitesimal, 
a v ,a. c > are real parameters to be defined, while the 
extended fluid fields {p, Y,pi} and the volume force f 
are all assumed independent of s. Hence, with respect to 
e they scale 



p , V,pi,f ~ o(e°) 



(57) 



As a result, for suitably smooth fluid fields (i.e., in va- 
lidity of Axiom 1) and appropriate initial conditions for 
fi(r,t), it is expected that the first requirement actually 
implies in the whole set x / the condition of closeness 
/j(r, t) = fl q (r, t) [1 + o(e)\ , consistent with the LB As- 
sumption #4. To display meaningful comparisons with 
previous LBM's let us introduce the further assumption 
that the fluid viscosity is small in the sense 



Ai~o(e a "), 



(58) 



with > 1 another real parameter to be defined. 
Asymptotic approximations for the corresponding LB- 
IKE [Eq. (l30|) ] can be directly recovered by introducing 
appropriate asymptotic orderings for^the contributions 
appearing in the source term S% = Si. Direct inspec- 
tion shows that these are provided by the (dimensional) 



M eff 



1 dp 

K f J = ~|v-n-vp|, 



p.b 

M^ ff = - | M V 2 V 



(59) 
(60) 
(61) 



The first two M**J and M* f / are here denoted respec- 
tively as (first and second) pressure effective Mach num- 
bers, driven respectively by the pressure time-derivative 
and by the divergence of the pressure anisotropy II — pi. 

Furthermore, My is denoted as velocity effective Mach 
number. Physically relevant examples [of asymptotic 
LBM's] can be achieved by introducing suitable orderings 
in terms of the single infinitesimal e for the parameters 
Mpff , M^/J and My^. We stress that these orderings, in 
principle, can be introduced without actually introducing 
restrictions on the fluid fields, i.e., retaining the assump- 
tion that the extended fluid fields are independent of e. 
Interesting cases are provided by the asymptotic order- 
ings indicated below. 



7A - Small effective Mach numbers (MpfJ,M^ f b f and 



A/ V J/ ) 



An important aspect of LB theory is the possibility 
of constructing asymptotic LBM's with prescribed accu- 
racy with respect to the infinitesimal parameter e, in the 
sense that the fluid equations are satisfied at least cor- 
rect up to terms of order o(s n ) included, with n = 1 
or 2, namely ignoring error terms of order o(e n+1 ) or 
higher. Let us, first, consider the case in which all pa- 
rameters M^/ , M^J and M^J are all infinitesimal w.r. 
to e (low- effective- Mach numbers). Since the parameters 
c and v c are free, they can be defined so that that there 
results c ~ v c ~ \Jo(e) [which implies ct c = a v = 1]. 
This requires 



M eff 



o(s 2 ). 



(62) 



If, we consider a low- viscosity fluid for which the kine- 
matic viscosity v = /x/p Q can be assumed of order e [and 
hence a M = 1] it follows that 



eff 



o(e 2 ). 



(63) 



Thanks to the assumptions (l54 |) -([58 |) there follows V • 
II — Vp ~ o(e) and /iV 2 V ~ o(e), which implies that 

the source term ignoring corrections of order o(e 2 ), 
becomes 



Si 

SAi 



Sm [1 

72 



o( £ )] 
ai • f . 



(64) 
(65) 
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It is immediate to determine the corresponding moment 
equations, which read: 



1 dpi 
c 2 dt 

NV = "<~ 2 



v-v = 

f o(e 2 ), 



0. 



(66) 
(67) 



Formally the first equation can be interpreted as an evo- 
lution equation for the kinetic pressure p\. Nevertheless, 
in view of the ordering ([62]) it actually implies the iso- 
choricity condition 



V-V = + o(e 2 ). 



(68) 



Instead, the second one [Eq. (f6"Tj) ]. due to the asymp- 
totic approximation (|63p . reduces to the Euler equation. 
Therefore in this case the asymptotic approximation (|64[) 
is not adequate. To recover the correct Navier-Stokes 
equation a more accurate approximation is needed, real- 
ized requiring that the hydrodynamic equations are sat- 
isfied correct to order o(e 3 ). A fist possibility is to con- 
sider a more accurate approximation for the source term. 
Restoring the pressure and viscous source terms in (|64[) 
there results the asymptotic source term 



01 

,.2 



dpi 
dt 



a, • (f!-^V 2 V) 



(69) 



where in validity of the previous orderings 

% = S Bl [1 + o(e)] . (70) 
The corresponding moment equations become therefore 



V- v = o, 

NV = + o(e 3 ). 



(71) 
(72) 



It is remarkable that in this case the isochoricity condi- 
tion is exactly fulfilled, even if the source term is not the 
exact one. For the sake of reference, it is interesting to 
mention another possible small-Mach-number ordering. 
This is obtained imposing for the parameters c and v c 



1 

i 



(73) 
(74) 



while requiring for v = fJ,/p a the same constraint adopted 
by asymptotic LBM's, namely Eq.([T7|). In this case one 
can show that the moment equation (|72|) is actually satis- 
fied correct to order o(e 3 ), while the isochoricity condition 
is only satisfied to order o(e 2 ). The following theorem 
can, in fact, be proven: 



Theorem 3 - Low effective-Mach-numbers asymptotic 
approximation 

In validity of THM.l, let us invoke the following as- 
sumptions: 



1) LB assumptions #3 and #4 f or the discrete kinetic 
distributions fi (i = 0, 8); 

2) the free parameters c and v c are assumed to satisfy 
the asymptotic orderings J7<S| j, f74P ; 

3) the fluid viscosity fi is assumed of order /i ~ o(e) 

4 ) the fluid viscosity fi is prescribed so that the kine- 
matic viscosity v — p/p is defined in accordance to 

5) the kinetic pressure p\ is assumed slowly varying in 
the sense 



dlnpi 
dt 



o(e). 



(75) 



It follows that the source term is approximated by 
Eg. \64^ and moment equations are provided by the 
asymptotic eguations: 



J_dpi 

c 2 dt 

NV = 1 -'- 3 



V- v = o 

f o(e 3 



o(e 3 ), 



(76) 
(77) 



i.e., the isochoricity and NS eguation are recovered re- 
spectively correct to order o(e 2 ) and o(e 3 ). 
Proof 

First we notice that the ordering assumptions 2)-5) 
require 



M eff 



o(e 2 ), 
( £ 4 ), 



(78) 
(79) 
(80) 



which imply at least the validity of Eqs. ([6~4]) -(f6"7l). The 
proof of Eqs. ([76|) and (|77|) is immediate. In both cases 
it sufficient to notice that in validity of hypotheses l)-3) 
and in terms of a Chapman-Enskog perturbative solution 
of Eq. (|30| there results actually 



-/iV 2 V - V- n + Vp = + o(e 3 ), 



(81) 



and hence Si reduces to Eq. ([64|) . 

The predictions of THM.3 are relevant for comparisons 
and to provide asymptotic accuracyestimates for previ- 
ous asymptotic LBM's [see Refs. [l7J, |2lj, |40j]. In fact, 
the asymptotic moment equations (|76[) and (|77[) formally 
coincide with the analogous moment equations predicted 
by such theories, when the kinetic pressure p is replaced 
by the fluid pressure p\ (i.e^ if the function P Q (t) is set 



identically equal to zero). |17j, l21|, |40( ■ Nevertheless, the 



accuracy of customary LBM's depends on the properties 
of the solutions of INSE. In fact, if one assumes 



91npi 
dt 



o(e°) 



(82) 



the customary (V,p) asymptotic LBM [17J, |2l|, |40j result 
actually accurate only to order o(e 2 ). Therefore, in such 
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case to reach an accuracy of order o(e 3 ) the approxima- 
tion (|69|) must be invoked for the source term. 

The other interesting feature of Eas. ([75|) and (|77|) is 
that they provide a connection with the artificial com- 
pressibility method (ACM) postulated by Chorin [o3j |. 
previously motivated merely on the grounds of an asymp- 
totic LBM [2lj. In fact, these coincides with the Chorin's 
pressure relaxation equation where c can be interpreted 
as sound speed of the fluid. However - in a sense - this 
analogy is purely formal and is only due to the neglect of 
the first pressure source term in Si. It disappears alto- 
gether in Eq. ([7Tj) if we adopt the more accurate asymp- 
totic source term (|69|) . A further difference is provided 
by the adoption of the kinetic pressure p\ which replaces 
the fluid pressure p (used in Chorin approach). We stress 
that the choice of p\ here adopted, with P {t) determined 
by the entropic principle, represents an important differ- 
ence, since it permits to satisfy everywhere in ft x / the 
condition of strict positivity for the discrete kinetic dis- 
tribution functions. 

7B - Finite pressure-Mach number MpfJ 

Another possible asymptotic ordering, usually not per- 
mitted by customary asymptotic LBM's, is the one in 
which the test particle velocity is finite, namely c ~ o(e°), 
the viscosity remains arbitrary and is taken of order 
/i ~ o(e ) while again v c is assumed v c ~ l/o(e 2 ) [i.e., 
a c = v c = 0, a u = 2]. In this case the pressure Mach 
M^fJ number results finite, while velocity and the sec- 
ond pressure Mach numbers are considered infinitesimal, 
respectively of first and second order in e, namely 

M;!J ~ o(e°), 

M$" ~ o(e), (83) 
M e p i f ~ o(e 2 ). 

To obtain the fluid equation with the prescribed accu- 
racy, say of order o(e 2 ), it is sufficient to approximate 
the source term Si in terms of Si = S^j [l + o(e 2 )] . The 
set of asymptotic moment equations coincide therefore 
with Eqs. (j71"j) . ([72"f . Again, the isochoricity condition is 
exactly fulfilled, while in this case the NS equation is 
accurate only to order o(e 2 ). 

TC - Small effective pressure-Mach numbers 
(M^J,M^ b f ) and finite velocity-Mach number (M^ ff ) 

Finally, another interesting case is the one in which the 
fluid viscosity /x remains finite (strongly viscous fluid), 
i.e., in the sense fi ~ o(e°) [i.e., a p = 0] while both 
parameters c and v c are suitably large, and respectively 
scale as c ~ l/o(e), v c ~ l/o(e 2 ) [i.e., a c = 1, a v = 2]. 
Due to assumptions one obtains V • II — Vp ~ 



o(e 2 ) and ^iV 2 V ~ o(e°). It follows that the effective 
Mach numbers scale respectively as 

M$ - o(s 3 ) (84) 

If we impose on \i also the same constraint set by Eq. (fTT)) . 
the customary asymptotic LBM's can be invoked also in 
this case. However, since the first pressure and veloc- 
ity Mach numbers are only second order accurate, the 
NS equation is recovered to order o(e 2 ) only. Never- 
theless, it is possible to recover with prescribed accuracy 
the fluid equations (jTT j) . (|72|) . This is obtained adopt- 
ing the source term Si = Ssi [see Eq. (|69|) ]. As a basic 
consequence, the isochoricity equation is satisfied exactly 
(hence no meaningful analogy with Chorin's approach 
arises), while the NS equation results correct to order 
o(e 3 ). These results provide a meaningful extension of 
the customary asymptotic LBM's. We stress that the 
entropic approach here developed holds independently of 
the asymptotic orderings here considered [for the param- 
eters MpfJ , ,M^J^]. Thus it can be used in all cases 
to assure the strict positivity of the discrete distribution 
function. 

8 - CONCLUSIONS 

In this paper we have presented the theoretical foun- 
dations of a new phase-space model for incompressible 
isothermal fluids, based on a generalization of customary 
lattice Boltzmann approaches. We have shown that many 
of the limitations of traditional (asymptotic) LBM's can 
be overcome. As a main result, we have proven that 
the LB-IKT can be developed in such a way that it 
furnishes exact Navier-Stokes and Poisson solvers, i.e., 
it is - in a proper sense - an inverse kinetic theory for 
INSE. The theory exhibits several features, in particular 
we have proven that the integral LB-IKT (see Sec. 5): 

1. determines uniquely the fluid pressure p(r,t) via 
the discrete kinetic distribution function without 
solving explicitly (i.e., numerically) the Poisson 
equation for the fluid pressure. Although analo- 
gous to traditional LBM's, this is interesting since 
it is achieved without introducing compressibility 
and/or thermal effects. In particular the present 
theory does not rely on a state equation for the 
fluid pressure. 

2. is complete, namely all fluid fields are expressed as 
momenta of the distribution function and all hy- 
drodynamic equations are identified with suitable 
moment equations of the LB inverse kinetic equa- 
tion. 
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3. allows arbitrary initial and boundary conditions for 
the fluid fields. 

4. is self-consistent: the kinetic theory holds for ar- 
bitrary suitably smooth initial conditions for the 
kinetic distribution function. In other words, the 
initial kinetic distribution function must remain ar- 
bitrary even if a suitable set of its momenta are 
prescribed at the initial time. 

5. the associated the kinetic and equilibrium distri- 
bution functions can always be chosen to belong to 
the class of non-Galilei-invariant distributions. In 
particular the equilibrium kinetic distribution can 
always be identified with a polynomial of second 
degree in the velocity. 

6. is non- asymptotic, i.e., unlike traditional LBM's it 
does not depend on any small parameter, in partic- 
ular it holds for finite Mach numbers. 

7. fulfills an entropic principle, based on a constant-H 
theorem. This theorem assures, at the same time, 
the strict positivity of the discrete kinetic distri- 
bution function and the maximization of the as- 
sociated Gibbs-Shannon entropy in a properly de- 
fined functional class. Remarkably the constant H- 
theorem is fulfilled for arbitrary (strictly positive) 
kinetic equilibria. This includes also the case of 
polynomial kinetic equilibria. 

A further remarkable aspect of the theory concerns the 
choice of the kinetic boundary conditions to be satisfied 
by the distribution function (Axiom II) and obtained by 
prescribing the form of the incoming- velocity distribution 
[see Eq.([36]>]. Thanks to Eqs.{34]),([35]), this requirement 
[of the LB-IKT] the boundary conditions for the fluid 
fields are satisfied exactly while the fluid equations are 
by construction identically fulfilled also arbitrarily close 
to the boundary. This result, in a proper sense, applies 
only to Dirichlet boundary conditions for the fluid fields 
[see Eqs.([H])]. Nevertheless the same approach can be 
in principle extended to the case of mixed or Neumann 
boundary conditions for the fluid fields. 

Moreover, we have shown that a useful implication of 
the theory is provided by the possibility of constructing 
asymptotic approximations to the inverse kinetic equa- 
tion. This permits to develop a new class of asymptotic 
LBM's which satisfy INSE with prescribed accuracy, to 
obtain useful comparisons with previous CFD methods 
(Chorin's ACM) and to achieve accuracy estimates for 
customary asymptotic LBM's. The main results of the 
paper are represented by THM's 1-3, which refer respec- 
tively to the construction of the integral LB-IKT, to the 
entropic principle and to construction of the low effective- 
Mach-numbers asymptotic approximations. For the sake 
of reference, also another type of LB-IKT, which admits 



as exact particular solution the polynomial kinetic equi- 
librium, has been pointed out (THM.lbis). 

The construction of a discrete inverse kinetic theory of 
this type for the incompressible Navier-Stokes equations 
represents an exciting development for the phase-space 
description of fluid dynamics, providing a new starting 
point for theoretical and numerical investigations based 
on LB theory. In our view, the route to more accurate, 
higher-order LBM's, here pointed out, will be important 
in order to achieve substantial improvements in the effi- 
ciency of LBM's in the near future. 



APPENDIX A 

The basic argument regarding the accuracy of the 
boundary conditions ado pted by customary asymptotic 
LBM's is provided by Ref. [4q . In fact, let us assume that 
on the boundary <5f£ the incoming distribution function 
fj \r w ,t) is prescribed according to Eqs.([33 | . ([37|) and 
([38]) . being ft \r w ,t) prescribed suitably smooth func- 
tions which are non vanishing only only for incoming dis- 
crete velocities for which (a,; — V m ) ■ n(r w , t) < 0. For 
definiteness, let us assume that f^~\Y w ,t) = f^ q (r w ,t) 
where f^ q (r w , t) denotes a suitable equilibrium distribu- 
tion. It follows that suitably close to the boundary the 
kinetic distribution differs from the Chapman- Enskog so- 
lution (|25jl . The numerical error can be overcome only 
discarding the first few spatial grid (close to the bound- 
ary) in the numerical simulation [46| . 



APPENDIX B 

Unlike standard kinetic theory, the distinctive feature 
of LB-IKT's is the possibility of adopting a non-Galilei 
invariant kinetic distribution function (i.e., non-invariant 
with respect to velocity translations). Here we re- 
port another example of discrete inverse kinetic theory 
of this type. Let us modify Axiom IV so that to permit 
that a particular solution of LB-IKE [Eq. ([30]) ] is pro- 
vided by fi = fl q . Here we identify f^ q with the (non- 
Galilei invariant) polynomial kinetic distribution defined 
by Eq. ([27|) but with the kinetic pressure p\ that replace 
the fluid pressure p. In this case one can prove that the 
source term Si reads 

Si = Sf 3 =Si + ASi, (85) 
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where 



AS, 



Wip 



1 



V a, 



V) • VV-- ai V-V 



3a, • V 



(86) 



a, V 



V 2 



Here Aq = N—p J^, where N is the Navier-Stokes opera- 
tor ([6|) , namely N\ is the nonlinear operator which acting 
on V yields AqV = p V-VV+V [ Pl - * (r)]+f!-/A7 2 V. 
Hence, invoking INSE, AS^ can also be written in the 
equivalent form 



AS,; = 



Wipo 



(a,-V)- VV--a,V-V 



w l d 



V -an 



3a, : • V 



(87) 



--WiPo&i ■ V 



a, V 



r.2 



The following result holds: 



Theorem Ibis - Differential LB-IKT 



In validity of axioms I-IV and the assumption that 
fi = fi 9 * s a particular solution of Eq. &30\) , the following 
statements hold: 



A) 1 



particular solution of LB-IKE [Eg. ftJOp J 



if and only if the extended fluid fields {V,pi} are strong 
solutions of INSE of class V29\), with initial and boundary 
conditions |0)-f5|), and arbitrary pseudo pressure Po (t) of 
class CW(J). 

Moreover, for an arbitrary particular solution fi and 
for arbitrary extended fluid fields: 

For an arbitrary particular solution fi : 

B) fi is a solution of LB-IKE [Eg. \30\) 1 if and only 
if the extended fluid fields {V,p\} are arbitrary strong 
solutions of INSE of class I129\), with initial and boundary 
conditions |0)-f3|), and arbitrary pseudo pressure Po (t) of 
class CW(7); 

C) the moment equations of L-B IKE coincide identi- 
cally with INSE in the set Q x I; 

D) the initial conditions and the (Dirichlet) boundary 
conditions for the fluid fields are satisfied identically; 

E) the source term Si is uniquely defined by 
Eqs.fSB.fSW; 

Proof: 

The proof of propositions A,B, C and D is analogous 
to that provided in THM.l. Assuming S, = S t (1) , the 
proof of B follows from straightforward algebra. In fact, 
letting fi(r,t) = f- q (r,t) for all (r,t) e x I in the 
LB-IKE [Eq.®], one finds that Eq.® is fulfilled iff 



the fluid fields satisfy the Navier-Stokes, isochoricity and 
incompressibility equations (P),© and ([3]). The proof 
of proposition E can be reached in a similar way. The 
uniqueness of the source term Sj is an immediate conse- 
quence of the uniqueness of the solutions for INSE. 
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